Matrix product state approach for a two-lead, multi-level Anderson impurity model 
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We exploit the common mathematical structure of the mimerical renormalization group and the 
density matrix renormalization group, namely, matrix product states, to implement an efficient 
numerical treatment of a two-lead, multi-level Anderson impurity model. By adopting a star-like 
geometry, where each species (spin and lead) of conduction electrons is described by its own Wilson 
chain, instead of using a single Wilson chain for all species together, we achieve a very significant 
reduction in the numerical resources required to obtain reliable results. We illustrate the power of 
this approach by calculating ground state properties of a four-level quantum dot coupled to two 
leads. The success of this proof-of-principle calculation suggests that the star geometry constitutes 
a promising strategy for future calculations the ground state properties of multi-band, multi-level 
quantum impurity models. Moreover, we show that it is possible to find an "optimal" chain basis, 
obtained via a unitary transformation (acting only on the index distinguishing different Wilson 
chains), in which degrees of freedom on different Wilson chains become effectively decoupled from 
each other further out on the Wilson chains. This basis turns out to also diagonalize the model's 
chain-to-chain scattering matrix. We demonstrate this for a spinless two-lead model, presenting 
DMRG-results for the mutual information between two sites located far apart on different Wilson 
chains, and NRG results with respect to the scattering matrix. 

PACS numbers: 78.20.Bh, 02.70.-c, 72.15.Qm, 75.20.Hr 



I. INTRODUCTION 

A very successful method for solving quantum impu- 
rity models is Wilson's numerical renormalization group 
(NRG)i"— . Recently, it has been pointed outii that the 
approximate eigenstates of the Hamiltonian produced 
by NRG have the structure of matrix product states 
(MPSs)<^ This observation established a structural re- 
lation between NRG and the density matrix renormal- 
ization group (DMRG)^"— because the states produced 
by the latter likewise have the form of MPS<^"— 

This structural relation between NRG and DMRG has 
opened up very interesting perspectives for combining 
advantageous features of both methods. In particular, 
the fact that DMRG, in essence, is a method for vari- 
ationally optimizing MPSs^ '^^i^^ can be used to devise 
a corresponding variational treatment of quantum impu- 
rity modelsi^ii^ This has the advantage that MPSs with 
much richer more complex structures can be adopted 
than those produced by standard NRG, entailing a much 
more efficient use of numerical resources. Concretely, the 
dimension D of the matrices from which the MPS is con- 
structed can be reduced very significantly, typically by 
several orders of magnitude. As a result, it becomes fea- 
sible to study complex quantum impurity problems that 
would be very challenging for standard NRG. 

In this paper, we illustrate this idea by calculating 
ground state properties of a multi-level quantum dot cou- 
pled to two spinful leads. Standard NRG treats the lat- 
ter as a single quantum chain with 2^ states per site (to 
account for two spin and two lead degrees of freedom), 
for which one typically needs D > 4000 to achieve sat- 
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Fig. 1: Quantum dot coupled to two leads. 



isfactory results. In contrast to the latter "single-chain 
geometry," we adopt here a MPS with a "star geometry," 
involving four separate chains, each with only two states 
per site, and variationally optimize one chain after the 
other. This enables us to obtain good results using ma- 
trices with D ranging between 16 and 36. This reduction 
in numerical memory resources relative to standard NRG 
illustrates the increased numerical efficiency alluded to 
above. Furthermore, we show that a numerically opti- 
mal basis, involving rotated Wilson chains, can be found 
by requiring that this new representation minimizes the 
mutual information between different chains. This opti- 
mal basis has an instructive physical interpretation: it is 
the basis in which the chain-to-chain scattering matrix is 
diagonal. 

This paper is structured as follows. In Sec.[n]we briefly 
review why standard NRG produces MPSs with a single- 
chain geometry and advocate the adoption of MPSs with 
an alternative star geometry. In Sec. IlIU we describe 
how a star-MPS representation of the ground state can 
be determined by variationally minimizing its energy. In 
Sec. lIVI we present proof-of-principle calculations of some 
ground-state properties and comparisons thereof to NRG 
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results. Finally, Sec. fVl illustrates how a numerically op- 
timal basis for the chains can be obtained by effectively 
minimizing the mutual information between two sites of 
different chains. 



II. MATRIX PRODUCT STATE ANSATZ 
A. Model 

We study a multi-level, two-lead Anderson impurity 
model described by the following Hamiltonian: 

H = H^ot + Hint + ^^Icads + -f^coupling, (1) 

where i^dot describes the eigenenergies of the m dot levels 

m 

i^dot = 51 XI ^^'^4sd^s: (2) 
Hint is the Coulomb interaction on the dot 

^int = | E 4sd.s4s'djs', (3) 

-fficads is the free lead Hamiltonian for A^i leads {a = 
l,...,iVi) 

ifleads = E^4o/feo.' 
kas 

and i?coupiing is the coupling between the dot levels and 
the leads 

^coupling = J2 + 4os^-) • 

ikas 

At a late stage of this work we became aware of work 
of Kashcheyevs et al suggesting to perform a singular 
value decomposition on -ffcoupiing which has the merit of 
decoupling some levels from some leads. Applying this 
idea to our system should also give some improvement in 
numerical efficiency. In general, however, all the levels 
will remain to be coupled to all leads. As we will show 
later, a more general scheme than just a singular value 
decomposition is capable of generating a new basis for 
the leads that will minimize the coupling of the leads 
amongst themselves. 

Following Wilson^ we adopt a logarithmic discretiza- 
tion of the conduction bands and tridiagonalize i/ioads + 
^^coupiing- As a result, the dot, represented by the "dot 
site," is coupled to the first sites of 2A'i separate "Wilson" 
chains, labeled by (a, s) 

^coupling ^WJ2 {fLsd^s + 4./0.s) (6) 

ias 

i^lcads = VK^i(l-KA-l) 

L-l"' , (7) 

X ^ A-t^„ (/L. /(n+l)as + H.C.) . 
n=0 



Here = (1 - A-"-l)(l - A^^n-l^-l^^ _ ^-2n-3y^ 

are coefficients of order 1, Tia — T^pV-^ the hybridization, 
p is the density of states, and 2W is the bandwidth of 
the conduction bands of the leads centered at the Fermi 
edge. We set the NRG discretization parameter A = 2 
throughout this paper. The length L of the Wilson chain 
determines the energy resolution with which the lowest- 
lying eigenstates of the chain are resolved. We typically 
choose L = 60. 

A standard NRG treatment of this model would com- 
bine all four Wilson chains into a single one, whose 
sites are labeled by a single site index k = 0, . . . , L [see 
Fig. [3ja)]. Each site would represent a 2^''^i-dimensional 
local state space, consisting of the set of states {|crfc)}, 
where the state label Cfc takes on 2^^' different val- 
ues. Then one proceeds to diagonalize the Hamilto- 
nian iteratively, as follows: suppose a short Wilson 
chain up to and including site fc — 1 has been diag- 
onalized exactly, yielding a set of eigenstates \ik) S 
span{{|cri)} (g) {|cr2)} (g) . . . S?) {|crfe_i)}}. Then one adds 
the next site, fc, to the chain, thereby enlarging the 
Hilbert space by a factor of 2^^' , diagonalizes the Hamil- 
tonian in this enlarged space, and truncates by discarding 
all but the lowest D eigenstates of the Hamiltonian. The 
latter can in general be written as linear combinations of 




Fig. 2; Iterative generation of matrix product states for a 
chain. 

the following form (illustrated in Fig. [2]): 

iz.+o = E4:i.N.)i-.)- (8) 

Iterating this procedure up to and including site L pro- 
duces eigenstates of the form 

NL+i)=4:i+,---4:i+j^/.)k^)-->L): (9) 

where sums over repeated indices are implied. Since such 
states are completely characterized by sums over prod- 
ucts of matrices, they have come to be known as matrix 
product states. The form of these MPS produced by 
NRG is analogous to the state for a chain as shown in 

Fig.m 

B. Star geometry 

One limiting factor for the accuracy of the NRG ap- 
proach is that a certain amount of information is lost at 
each iteration step due to truncation. In general, for a 
system with Ni bands (in the two-lead case which we will 
investigate below, Ni — 2), the dimension of the effective 
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Hilbcrt space is enlarged from D to D2^^' upon adding 
a new site to the Wilson chain. Thus, the larger Ni, the 
more information is lost during the subsequent trunca- 
tion of the Hilbert space back to dimension D, and the 
less accurate the NRG treatment is expected to be. 



LSS size: 2'^"' 
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Fig. 3: (a) Single chain geometry: a single Wilson chain of 
local dimension 2^^' coupled to one dot site of local dimension 
2^™ . (b) Star geometry: 2Ni Wilson chains (here iVi = 2 and 
a = l,r) , each with local dimension 2, coupled to two dot 
sites of local dimension 2™ . 

The main goal of the present paper is to illustrate that 
a very significant improvement of efficiency can be ob- 
tained as follows: instead of combining all 2N\ chains into 
a single Wilson chain of local dimension 2^^' ("single- 
chain geometry"), we shall treat them as separate chains, 
each with local dimension 2 and each coupled to the same 
set of dot levels ["star geometry," see Fig. [3jb)]. Al- 
though the total number of sites thereby increases from 
0{L) to 0{N\L), the dimension of the local state space 
per site is reduced from 2^^' to 2. We find that, due to 
the latter fact, the dimension D of the constituent ma- 
trices in the star-MPS can be chosen to be significantly 
smaller than in the chain-MPS. 

The change from single-chain to star geometry, how- 
ever, necessitates a change in truncation strategy for the 
following reason: in contrast to the single-chain geome- 
try, where each site represents a definite energy scale, in 
the star geometry a given scale is represented by a set 
of 2N\ sites, one on each of the star's chains, i.e., at lo- 
cations that are widely "separated" from each other on 
the star. Therefore, a truncation scheme based on energy 
scale separation, such as that used by standard NRG, can 
no longer be applied. Instead, we shall simply minimize^ 
the expectation value of the Hamiltonian within the space 
of all MPSs with the same star structure. This can be 
done efficiently by optimizing the matrices in the star- 
MPS one site at a time, and sweeping through all sites 
until convergence. 

To be explicit, we construct our star-MPS for the two- 
lead system as follows. In total 4 ~ 2N\ {N\ = 2) Wilson 
chains are connected to the dot. Each of these chains is 
very similar to the NRG MPS from above, except that the 
local state space (LSS) is only of dimension 2. To simplify 
the notation we drop the labels a and s whenever possible 
and incorporate them into the site index k, which from 
now on will be taken to uniquely determine a site in the 
whole star structure, still labels the LSS at site k. 



With this every Wilson chain can be represented as (see 
Fig. HI) 

|oo) ^ 4?o\4?ol • ■ • Wi)W2) ■ . ■ \<tl), (10a) 

where |(?) = |cri)|cr2) ■ • • I^l)- Here the label o stands for 
"outer," for reasons that will become clear below. We 
introduce an intuitive graphical representation for these 
MPS. Every A will be represented by a box and every 
index of A is depicted by a line attached to the box. For 
matrix products or other index summations the corre- 
sponding lines are connected. Using this representation, 
a single chain can be depicted as in Fig. [4l 



A, 



- A. -...-A 



- A, 



Fig. 4: Graphical representation of Eq, 



The fact that the Hamiltonian does not contain terms 
that flip spin up to down or vice versa suggests represent- 
ing the dot state space by two separate sites, representing 
all dot states having spin up or down, respectively [see 
Fig-Elb)]. Correspondingly, we also introduce two types 
of dot matrices, AI'^otI and A^'^«^\ which carry an extra 
index v that is being summed over to link the spin up and 
down subsystems. So we arrive at the starlike structure 
of Fig. [5] with two linked dot matrices (one for each spin) 
and two leads (left and right) attached to each: 



Oil 



(11) 



on 



This starlikc structure basically consists of two y- 
junctions, as discussed by Guo and White^^ next to each 
other. 

Hiding the explicit structure [Eq. (fTT|) ] of the MPS as 
illustrated in Fig. [51 we can write a state symbolically as 



(12) 



We call Eq. (fT^ the global representation of 

An important point to note is that this system is still 
effectively one dimensional, in the sense that if we cut out 
a given site, the system breaks apart into two (or three 
in case of a dot site) disjoint parts. We shall call the one 
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Fig. 5: MPS representation for a quantum dot coupled to two spinful leads. The lead chains are combined to big boxes for 
clarity. The indices of the dot matrices are labeled explicitly. 



containing the dot sites the "inner" part, the other one 
the outer part. As a consequence, it is possible to also 
give a "local" description of of the form 

\i,)^A^±\^,)\a,)\o,), (13) 

where {|crfe)} represents the LSS of the chosen site, {\ik)} 
is an orthonormal set of states representing the inner 
state space (ISS), namely, the inner part of the star with 
respect to the chosen site fc, and {|ofc)} is an orthonormal 
set of states representing the outer state space (OSS), 
namely, the outer part of the star. 

III. VARIATIONAL SITE OPTIMIZATION 
SCHEME 

We will use the MPS of Fig. [5] as an ansatz for the 
ground state of our system. In order to find the ground 
state we need to calculate the MPS lip) that minimizes 
the energy E = {ijj\H\ip) with the constraint of keeping 
the norm of {"ip) constanti^ Using A as Langrange mul- 
tiplier ensuring normalization we arrive at the following 
minimization problem: 

min((V'|i?|V') - A(?MV')) • (14) 

The key idea of the variational MPS optimization is to 
optimize every single A- matrix of |^) separately until the 
ground-state energy has converged. Therefore we insert 
the local MPS description from Eq. ^ into Eq. ^ 
and obtain 

mm (A51*i?(..^o').(,:.,o)4?l - a4?1*4?1) , (15) 

where -ff(i'o'crj,).(jo(TA,.) ^re the Hamilton matrix elements 
in the current effective bases 

i^(..'<o'),(«^.o) = W\{al\{z'\H\i)\ak)\o). (16) 

By setting the derivative of Eq. (flSl) with respect to the 
matrix elements of A^. to zero and replacing Xhy Eq, we 
obtain the following eigenvalue equation for Af^: 

H{i'a'^o'),iicrko)A'fo''^ = EqA^J^} . (17) 

The eigenvector with the smallest eigenvalue is the solu- 
tion to our minimization problem. So after having solved 



this eigenvalue problem for the current site k we replace 
Ak with the newly found eigenvector and move on to the 
next site in order to optimize that Ak' ■ We repeat the 
whole process (sweeping) until the ground-state energy 
has converged (see below). 

By following this procedure we succeed to divide a very 
high dimensional minimization problem into manageable 
smaller units. For general problems this can be a very 
bad approach as one can get stuck in a local minimum 
during the optimization. However, it has proven to work 
reliably when the site-site coupling varies smoothly and 
monotonously. In our case the Hamiltonian has only 
nearest-neighbor interactions and there are no long-range 
correlations in the system. As a result, the system reli- 
ably converges without getting stuck in local minima. 



A. Updating the A matrices and changing the 
eflFective basis states 

When updating A matrices during sweeping, one must 
ensure that two conditions are satisfied. First, whenever 
we use the local description of Eq. (|13p . we rely on the 
basis states being orthonormal: (o/jjoj.) = Sg^ o'^. This 
condition translates to 

^K'l^k.'lt = I for fc' > (18) 

for all outer matrices with respect to site k. We will 
focus here on the OSS basis, everything works completely 
analogously for the ISS basis. 

Second, we also want to create an effective basis that 
spans a DMRG optimal Hilbert space, i. e., the states we 
keep for an effective basis are to be the ones having the 
largest weights in the density matrix of the current state 
(as described below). 

For definiteness, we consider an inward sweep and fo- 
cus on how to move the "current site" from fc to fc — 1. 
We assume that a new set of A matrices for site k has 
been obtained by energy minimization. The question is 
how to ensure that both above mentioned conditions are 
satisfied. As all the inner A matrices of site k — 1 have not 
changed since we optimized site k ~ 1 the last time when 
moving outwards, we only need to create a new effective 
OSS basis \ok-i) for site fc — 1. 
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Starting from the density matrix in the local descrip- 
tion of site k, 

p^^^ = \m\ = (19) 

suppose one traces out the inner part of this system to 
obtain reduced density matrix of the outer part and site 
k. 



(20) 



which corresponds precisely to the outer part with re- 
spect to site fc — 1 . 

Now employ the singular value decomposition (SVD) 
A = U SV"^ which exists for every rectangular matrix 
A. iS* is a diagonal matrix containing the singular values 
ordered by magnitude; U and V'^ are column and row 
unitary matrices, respectively, and obey U^U ~ V^V ~ 
1. Combine \ak) and \ok) to \lk) — \crk)\ok) and insert 
the SVD for Au 



Uim Smj {V^ )jl 



(k) 
^rcd 



(21) 



The second line follows since is diagonal, and we 
wrote p'j''^ = and \jk) = Vji\lk)- We see that the 
SVD automatically diagonalizes the reduced density ma- 
trix with the states ordered according to their weight. 

So all we actually have to do for moving the actual 
site from fc to fc — 1 is to calculate the SVD of the newly 
optimized Ak = USV'' . We then replace Ak ^ Ak ^ 
and Ak-i — > A^^i — Ak-iUS as illustrated by Fig. [H] 
By doing so we do not change the total state, since the 



US 



42t 



Fig. 6: Procedure for moving the actual site from fc to fc — 1. 
The matrices that are not orthonormalized in any direction 
are printed with gray background. The gray lines within the 
boxes indicate whether the row or column vectors are or- 
thonormal (with the local level associated with row or column, 
respectively). 



The relation between the singular values and the 
weights of the reduced density matrix can be used to 
optimize our choice for the dimensions of the respective 
effective Hilbert spaces: instead of using the same di- 
mensions for all A matrices in the system, which turns 
out to be inefficient for inhomogeneous ones like ours, 
we adopt as truncation criterion the demand that the 
minimum value of 5*^ at a given site is to be smaller 
than some threshold uimin (in our case typically taken as 
10~^). After calculating the singular values, we choose 
the matrix dimensions Dk at the corresponding bond fc 
(between site fc and its neighbor in the direction of the 
dot) according to the following recipe. We choose Dk 
large enough to ensure that the minimal singular value 
Smin(fc) fulfills s,^jj;jj(fc) < Wmin, but subjcct to this Con- 
straint choose Dk to be as small as possible, in order to 
minimize computational resources. 

It is instructive to also explore the relation between Dk 
and the bond entropy Sk of site fc, which can be computed 

(k) 

from the reduced density matrix p]^.^l^ at site fc according 
to 



Sk 



(24) 



The entropy Sk is a measure for the entanglement be- 
tween the traced out part of the system and the part kept 
in the description of Thus, large Sk implies large 

Dk, which turns out to be roughly proportional to e^'' . 
The dimensions Dk resulting from the above criterion 
for the singular values Sniin(fc') together with the expo- 
nentiated bond entropy e^'' associated with the reduced 
density matrix at bond fc are shown in Fig. [71 This figure 
shows, first, that a larger dimension is required near the 
dot, and second, that e^'' (times a constant) is a rather 
good indicator of the required dimension Dk- For the 
limiting case of a reduced density matrix with uni- 
form weights Pj{k) = 'Wii^i ^ [Ij^fcli the exponentiated 
bond entropy then gives e"^*-' = Dk- Thus, Dk is a upper 
bound to e"^'= ^ The dip at fc = for the bond between 
the two spin subsystems (dimension Dy) is due to the 
fact that there is only a density-density interaction along 
this bond but no particle exchange. For our system we 
found that it is sufficient to have dimensions of 36 or less 
near the dot. 



product 

Ak-iAk = Ak-iAk (22) 

remains unchanged. Thus we create an effective or- 
thonormal OSS basis, 

\ok-,) ^ M^-l,^\ak)\ok) (23) 

which at the same time is DMRG optimal. 

The so-called site optimization procedure outlined 
above, where we optimize the A matrices directly, is 
equivalent to one-site finite-size DMRG. 



B. Sweeping sequence 

In principle the order in which we optimize the single 
matrices during a sweep is not important. However, it 
is both convenient and more efficient to move only to a 
neighboring site (and not further) for the next optimiza- 
tion step. In this way we need to change the actual site 
only by one in order to get the desired new local descrip- 
tion. Having our MPS ansatz structure in mind, this 
requirement immediately suggests a particular order of 
sweeping, shown in Fig. [51 Starting from the far end of 
any chain we move in toward the dot matrix and then out 



6 




* k 



Fig. 7: (Color online) The solid line shows the dimension 
needed at bond k of the spin up chain to satisfy Wmin = 10~® 
for the reduced density matrix at each bond (negative k cor- 
respond to the left chain). The dashed line displays the ex- 
ponentiated bond entropy 6^*= multiplied by 4.5 to visually 
match the -Dfc,min curve for large k. Here A: = corresponds 
to the "vertical" bond between the two spin subsystems. The 
two insets show spectra of reduced density matrices at dif- 
ferent bonds k indicated by the vertical dashed lines of the 
main plot. The data shown in this figure has been obtained 
from the ground state of the four-level model shown in Fig. [10] 
with e = —1.7U. In general, the maximum dimension needed 
depends strongly on the model parameters. 



again along another chain. We repeat this until wc have 
covered the whole system. Sweeping that way (solid blue 
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Fig. 8: (Color online) Sweeping sequence. For clarity we place 
the spin up and spin down parts on top of each other to 
emphasize the star-like structure. The solid blue line depicts 
the standard sweeping sequence. 

line in Fig. [5]) we optimize the two dot matrices three 
times but all the other sites only twice. If one wants to 
optimize all sites twice during a sweep one can once skip 
the optimization step at the dot sites as indicated by the 
dashed blue line. 



As the dot matrices are by far the biggest in the sys- 
tem, optimizing them takes much longer than optimizing 
any of the chain matrices. Thus by skipping the dot op- 
timization step once, we can reduce the computational 
time needed for a single sweep. However, since the dot 
optimization step also has the biggest effect for improv- 
ing our MPS ansatz, skipping its optimization once has to 
be compensated by performing more sweeps to achieve as 
good convergence of the ground state as in the case where 
we perform three optimizations at the dot matrices. We 
compared both approaches for our model and found no 
significant differences in the overall performance. 

We stop the sweeping when the MPS has converged. 
To probe the convergence we compare the MPSs before 
and after sweep N, IV'at-i), and \iPn)- If the change in 
overlap. 



1 - |(-0Ar-l |^Ar)| < £, 



(25) 



is smaller than a certain threshold, we stop the sweeping. 
We typically use e = 10 ~^ and need 10-15 sweeps. This 
depends crucially on the system parameters, though, and 
in some cases we need to perform up to 25 sweeps. 



C. Numerical costs 

The most computational effort is needed for solving 
the eigenvalue problem [Eq. (fT7|) ] for the minimal eigen- 
vector. We use the Lanczos method for solving Eq. (|17p . 
which is an iterative method and requires the calculation 
of Hltp) in the local picture once for every iteration. As 
we cannot influence the number of Lanczos iterations in 
our implementation, we will only investigate the costs 
of calculating H\ip), which are given by the costs of the 

matrix-matrix multiplication Z^jocr^ ^(^'o'<T^,),(^o<Tfc)^lo'°'■ 
The costs of a matrix-matrix multiplication is given by 
the size of the outcome times the dimension of the index 
being summed over. -ff(i'o'CT;.),(iocrfc) splits up into a sum 

of different terms, such as (cD^'^ak ® {ck+i)o'oi each con- 
sisting of a direct tensor product of operators living in the 
ISS, OSS or LSS. Thus the product H\ip) can be spht up 
into smaller matrix products. By looking at the struc- 
ture of the Hamiltonian ([1]), one recognizes that there 
will be no terms containing tensor products of operators 
from the ISS and OSS, since they would correspond to 
next-nearest-neighbor terms, but tensor products with 
one operator from the LSS and the other one from the 
ISS or OSS. These terms lead to multiplications over an 
index of length Dd, being the product of the dimensions 
of the ISS and LSS. If the current site is the dot site, the 
size of the resulting matrix is D^Dyd'" and thus the costs 
for a single multiplication H\ip) at a dot site is given by 



dot 



0{D^D^(r"Dd) = 0{D^Dyd"'+^). 



(26) 



In case of a chain site instead of a dot site exactly the 
same reasoning applies and because of the smaller matrix 
size the costs reduce to 0{D^d^). From Eq. ([^S]) we see 
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that optimizing the dot sites is the most expensive step 
in the optimization and scales particularly unfavorably 
when the number of dot levels m is increased. 



D. Bond optimization 

As an alternative to the site optimization scheme dis- 
cussed above, we can begin to move the current site as 
in Fig. El to obtain Ak-iiUS)V'', where Ak = USVl 
At this step we can represent the overall state as \^) — 
(f/S')ij.ofc_i |jfc)|ofc_i). Now we perform the optimization 
on B = US in complete analogy to the site optimiza- 
tion and obtain a new B. Then A^-i is replaced by 
Ak-i = Ak-iB which results in a state with the actual 
site k — 1. We call this process "bond optimization" as 
the matrix we actually optimize is somehow located at 
the bond between two original sites. 

One can easily see that the costs for calculating 
H{i'o').{io)Bio arc 0{D^) and thus independent of the 
number of dot levels. Considering only the costs for a 
single sweep the bond optimization scheme will be con- 
siderable faster than site optimization, which is especially 
expensive at the dot sites. This advantage, however, is 
compromised to some extent by the slower convergence 
of the bond optimization due to the optimization taking 
place within in a much smaller effective Hilbert space. 
This makes more sweeps necessary and also enforces a 
lower threshold in Eq. ([25)) as convergence criterion. It 
turned out to be very difficult to judge the convergence 
of the bond optimization scheme based on Eq. ([25]) es- 
pecially if one starts from a state not too different from 
the actual ground state because in such cases the conver- 
gence can be really slow and one might wrongly consider 
the state already converged. 

However, one might try to avoid unnecessary site op- 
timizations at the beginning of the sweeping and use 
cheap bond optimizations instead and switch after sev- 
eral sweeps to the site optimization scheme to make use 
of the better convergence properties. 

IV. RESULTS FOR LOCAL OCCUPATIONS 

We used the approach described above to calculate the 
ground state and level occupancies of a spinful multi- level 
quantum dot coupled to two leads. Throughout this part 
we fix the Coulomb interaction U ~ Q.2W, 2W being the 
bandwidth, and use the convention W = 1. 

The results shown below demonstrate that it is possi- 
ble to calculate local ground state quantities of a com- 
plex quantum dot efficiently using this approach. Al- 
ready with calculating the occupation of the dot levels it 
is possible to investigate the stability diagram of small 
quantum dots.— Under certain conditions, local occu- 
pancies can be related to phase shifts, which in turn can 
be used to calculate the conductance through a quantum 
dot^ 
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Fig. 9: (Color online) Dot level occupation for a spinless two- 
level system, with ei,2 = e ± A/2, level spacing A = 0.117 
and couplings Ti; = Ti^ = 0.005f/, Ta; = T2r = SOTu- This 
parameter set was used in Sindel et al. (Ref. [20l ) N = ^ {m + 
71,2) is half the total dot occupation. Note that the sign in Via 
just serves as an indication of the sign of the related hopping 
matrix element Via in the Hamiltonian. 



First we consider the simpler case of a spinless two- 
level model with level positions 61^2 = e ± A/2, coupled 
symmetrically to two leads. NRG works very reliable for 
this kind of impurity model. The lower of the two levels is 
assumed to couple significantly stronger to the leads. We 
calculated the occupation, Ui = {d\d^), of both levels as 
a function of e, using both our MPS approach and NRG. 
In Fig. [9] we show the occupation of both levels as we 
sweep the gate potential by shifting the levels from be- 
low towards the Fermi edge of the leads and then further 
above. At the beginning of this process mainly the lower 
level starts to empty. This is due to the much bigger cou- 
plings r2 of the lower level compared to the upper level 
and results in an occupation inversion situation where the 
energetically higher level has higher occupation than the 
lower level. A second consequence of the small couplings 
Fi is the sharp transition of the occupation of the upper 
level from almost filled to almost empty. Once the upper 
level is almost empty the dot system may gain energy 
by increasing the occupation of the lower level without 
having to pay Coulomb energy. This leads to the non- 
monotonic occupation of the lower level, known as charge 
oscillation. See Sindel et al.— for a more detailed discus- 
sion. The results for the level occupation of the simple 
spinless model as shown in Fig. [9l demonstrate excellent 
agreement between both NRG and DMRG calculations. 
The relative difference of the ground-state energies ob- 
tained by NRG and MPS was on average 10~^. 

We demonstrate the power of the MPS approach by 
considering a spinful four-level dot coupled asymmetri- 
cally to two leads, a system sufficiently complex that 
its treatment by NRG is a highly challenging task. We 
therefore have no NRG reference data for this system and 
present only DMRG results. For every dot level we cal- 
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Fig. 10: (Color online) Dot level occupation for a spinful 
four-level system. We parametrize the dot level energies as 
Eia = e + ± _B/2 for s =t, i, where B represents the ap- 
plied magnetic field with B = 0.2U and e a gate voltage, with 
= (-0.1,-0.03,0.07, 0.1)(7 The coupling of the dot levels 
are chosen asymmetrically Fi,. — SiVa with s; — (1, —1, —1, 1) 
and Til = (0.5,0.02, 1,0.7)0.217. 



culatc the occupation riis ~ {d-l^d^^) as a function of gate 
voltage, as shown in Fig. [TOl This calculation is solely 
performed within the site optimization scheme. We kept 
the effective dimensions for all A matrices describing the 
leads the same compared to the two-level plot, only the 
LSS size at the dot matrices was increased, thus demand- 
ing more computational time for the optimization at the 
dot. 

For the four-level system we chose random values for 
the level couplings F varying over two orders of magni- 
tude. Moreover, as the couplings have been chosen asym- 
metric, one cannot simplify the model by decoupling cer- 
tain linear combinations of the leads, while keeping the 
remaining relevant degrees of freedom. The occupation 
of the individual levels shows very rich behavior. By 
sweeping the gate potential similar to the spinless case 
above, we find the sharpest transition for the second level 
(ri2t, 7^24.)- The couplings of this level are one magnitude 
smaller than all other couplings causing this sharp tran- 
sition and associated with it charge oscillations in all the 
other levels. 



V. ROTATION TO OPTIMAL BASIS OF 
WILSON CHAINS 

As described above the use of a star-shaped MPS 
works well for local quantities. However, one might 
ask the question whether introducing such a geometry 
causes a loss of longer-ranged correlations between dif- 
ferent chains. To be able to assess this question we con- 
sider two sites in different chains c ^ c', both at distance 
k from the dot. The mutual informationii 1^'^ (k) con- 
tained between these two sites is given by 
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Fig. 11: (Color online) (a) Mutual information ij, 
two sites situated in different leads but at equal distances k 
from the dot, for a spinless four- level, two- lead model with dot 
levels ei/U = (-0.1, -0.03, 0.07, 0.1) + e, e = -2U fixed, cou- 
plings Tir = (0.3,-0.02,-1,0.2) and Ta = (0.5,0.08,1,0.7) 
and A = 3. The dashed line shows Ip'~ for the system with 
the leads in the original basis of Eq. ([1]), whereas the solid 
line shows after the leads have been rotated by the (fixed, 
k independent) optimal angle ^opt obtained from Fig. I12f a'). 
(b) Exponentiated bond entropy e®'' along the right chain 
of the system both prior (dashed line) and after (solid line) 
the rotation with Sopt, indicating an effective reduction in the 
required matrix dimension Dk close to the impurity for the 
rotated system by about i for the same numerical accuracy. 



with the entropy S 



Sp = -ti (plnp) . 



(28) 



5„ 



(27) 



Here p'j.^^{k) is the reduced one-site density matrix ob- 
tained by tracing out the entire system except for site 
k in chain c. Likewise p'ji^^{k) is the reduced two-site 
density matrix, obtained by tracing out all sites except 
two. situated at a distance k from the dot in two differ- 
ent chains, c and c'. Ip'^ (fc) is a measure for how much 
information the sites contain about each other. As a con- 
sequence, a decaying Ip'^ (k) as a function of distance k 
indicates that chains c and c' effectively decouple. 

For simplicity and to make a comparison with NRG 
feasible, we restrict ourselves to the spinless case, e.g.. 
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we only look at the spin-up part of the original four- 
level system, however with different couplings compared 
to the parameters used for Fig. [TUl As NRG treats both 
the left and right lead in a combined single chain we can, 
nevertheless, study the effect of "unfolding" the two parts 
of the NRG chain. 

If we calculate 7^''' for this spinless two-lead Hamilto- 
nian as it stands, the correlations between two sites on 
opposite sides of the dot but at equal distance from it 
are found to decay only very weekly with k [Fig. Ilir a). 
dot-dashed line]. This illustrates, on the one hand, that 
our MPS ansatz does successfully capture correlations 
between sites representing comparable energy scales, in 
spite of the fact that in the star geometry they lie "far" 
from each other (namely on different chains). On the 
other hand, it also raises the question whether one can 
choose a (numerically) better suited basis for the leads 
that effectively does decouple different chains far from 
the dot. Since in that case the correlations would intrin- 
sically decay with distance from the dot, less numerical 
resources would be required to capture all correlations 
accurately. 

Indeed, we shall show that it is possible to choose such 
an optimal basis by making a suitably chosen unitary 
transformation which rotates the lead degrees of freedom 
into each other in an "optimal" way to be described be- 
low. When the leads are first rotated by a certain optimal 
angle of rotation 6'opt (defined precisely below) and ij;^' 
is calculated in this rotated basis, then 7^''' is found to 
decay rapidly with fc, see solid line in Fig. [TlTa). 

We begin with the observation that the labeling of the 
unfolded chains with a = l,r is arbitrary. We can choose 
any linear combination of I and r as new basis, e.g., for 
symmetric couplings to the dot it is well known that with 
the symmetric and antisymmetric combination only the 
symmetric lead couples to the dot while the antisymmet- 
ric lead is completely decoupled. To be specific, we can 
introduce a unitary transformation acting on the original 
lead states specified in the Hamiltonian 



f/Sna- — Upafa 



(29) 



independent of the site n and spin ct, acting only on the 
lead index a. For systems with time-reversal symmetry, 
the unitary matrix is always chosen real. So in our case, 
for iVi = 2 [/ = U{9) is a real two-dimensional matrix 
and can be thought of as a planar rotation parametrized 
by a single angle 9. The optimal basis for DMRG treat- 
ment would have minimal correlations between the ro- 
tated chains. The angle of rotation 6 can be restricted to 
9 6 [0,7r/2] as we choose to ignore the particular order 
and relative sign of the new basis vectors. In order to find 
the optimal angle it is sufficient to look at the reduced 
two-site density matrix pl'^^{k). As the Hamiltonian dT]) 
preserves particle number, this density matrix is a 4 x 4 
matrix in block form: a 1 x 1 block for both the zero- 
particle and two-particle sectors and a 2 x 2 block for the 
one-particle sector. 
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Fig. 12: (Color online) Optimal basis for the leads of a spinless 
four-level, two-lead system (same parameters as for Fig. 1111 
but with varying e). (a) Optimal angle of rotation ^opt 
for the leads obtained by diagonalizing p^'^^(k = L) for the 
DMRG calculation (red symbols) in comparison with angle 
that diagonalizes the scattering matrix calculated with NRG 
(blue line), ^opt is defined mod 7r/2. (b) Dot level occu- 
pation. A'' = i rii is the rescaled total dot occupation. 
Rapid changes in the angle 9apt coincide with rapid shifting 
of dot-level occupations, (c) Truncation error (accumulated 
discarded density-matrix eigenvalues) of the DMRG calcula- 
tion considering two neighboring sites at a time for a rotated 
and non-rotated system. We typically used 20 sweeps for the 
DMRG calculations. The truncation error is significantly re- 
duced for the rotated system except for the points where ^opt 
actually shows rather rapid transitions through Sopt = it- 
self. At these points the leads are already decoupled from the 
outset. 



Finite off-diagonal elements of this 2x2 block show 
that both sites are effectively correlated with each other. 
However, by diagonalizing this block of ^'^^(fc) via a real 
unitary matrix U we immediately obtain a rotated lead 
basis according to Eq. ([29|) . So the angle of rotation 
^opt can be found by diagonalizing p^'^^{k). It is most 
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desirable to decouple the far ends of the chains best, so 
we choose 9 = 9{k = L), where U[6{k = L)] diagonalizes 

By applying the transformation U{9) to the Hamilto- 
nian ([1]) only the tunneling elements to and from the dot 
levels are changed 

V/3icr = U {9 opt) PaVaia- (30) 

This way, we have obtained a new lead basis for our 
Hamiltonian that is better suited for the DMRG calcu- 
lations, as long ranging correlations are suppressed in 
this basis. As we benefit already from a rotation in the 
leads even if the angle is only close (but not equal) to 
the optimal choice ^opt , it is feasible to start with a small 
system (of only, say, 14 sites per Wilson chain) in order 
to obtain an approximate value for ^opt; the latter can 
then be used to rotate the leads of a bigger system, from 
which a better determination of the optimal angle can be 
extracted. 

In Fig. [m we show the optimal angle of rotation 6'opt 
for a spinlcss four-level system. We compare with NRG 
calculations where we diagonalize the T-matrix 

r„^= lim (v^G''°\uj)v) , (31) 

where G'^°^ is the local retarded Green's-function matrix 
calculated by standard NRG techniques^i and V is the 
tunneling matrix from the Hamiltonian. The angle ex- 
tracted from the diagonalization of the T matrix [i.e., 
from requiring that U{9)TU^{9) be diagonal] is shown 
as a solid line in panel (a). Remarkably, this line agrees 
quantitatively with the 6'opt values found by DMRG. This 
shows that the angle of rotation that minimizes correla- 
tions between the two rotated leads has a clear physical 
interpretation: it also diagonalizes the scattering matrix, 
a result that is intuitively very reasonable. We note, 
though, that this fact cannot be used to determine ^opt 
before doing the DMRG calculation, as with the knowl- 
edge of the scattering matrix we would have already 
solved the system. Nevertheless, shorter systems can al- 
ready give a clean indication of the angle that decouples 
the chains. 

In Fig.lTTjwe demonstrate that by rotating the leads to 
the new optimal basis as suggested above it is possible, 
indeed, to ensure that lead degrees of freedom on differ- 
ent (rotated) Wilson chains become effectively decoupled 
from each other further out on the chains. Also the bond 
entropy Sk is reduced. If the leads are rotated into the 
optimal basis the mutual information drops quickly along 
the chains [see Fig.fT^a)]. and the truncation error is sig- 
nificantly smaller [see Fig. [T2l^c)], thus making numerical 
treatment less demanding. Note that rapid changes in 
the angle ^opt coincide with rapid shifting of dot-level 
occupations [see Fig. [T^ b)]. 



VI. SUMMARY 

Using the DMRG approach gives us the possibility to 
choose a more flexible MPS geometry compared to NRG. 
While in NRG one is bound to a simultaneous treatment 
of a single combined Wilson chain due to the require- 
ment of energy scale separation, this restriction can be 
lifted in a DMRG treatment. In our case of a two- lead 
Anderson model we modeled each spin and each lead by 
a Wilson chain on its own treated separately from each 
other. Thus we achieved a significant reduction in both 
the dimension of the LSS and the dimension D of the ISS 
and OSS at each site. The Hilbert space of one site in the 
single chain geometry is equivalent to the direct product 
of the Hilbert spaces of the 4 = 2N\ corresponding sites of 
each chain. So in order to map a star geometry descrip- 
tion into an equivalent single chain description, in the 
sense that the effective Hilbert spaces at every site have 
the same dimension, the dimension D' of the single chain 
A matrices would scale exponentially with the number of 
leads, D' ~ D^^\ as a consequence of the tensor prod- 
uct of the 2Ni smaller star geometry A matrices. Thus, 
adopting the star geometry reduces the numerical costs 
for treating the leads by D ~ j^'^/^^i^ Although this 
strategy has the consequence that the cost of treating the 
dot site increases significantly [see Eq. ([26)) ] , for all cases 
studied in this paper the latter effect is far outweighed 
by the decrease in costs for treating the leads. Indeed we 
found that dimensions D < 36 suffice for getting an accu- 
rate description of the system (note that when translated 
into a single chain this would result in a huge effective 
dimension of D' ~ D"* = 1.7 x 10^). 

Due to the fact that the Anderson Hamiltonian under 
consideration features only a density-density interaction 
term at the dot between electrons with different spin, the 
dot matrices can be conveniently split into two sets, one 
for each spin, yielding another gain in efficiency. As it 
turns out. the dimension Dy connecting to two sets of 
dot matrices can be chosen significantly smaller than D 
(cf. Fig. [71). In addition an optimal basis (in terms of 
numerical efficiency) for representing the leads has been 
determined, which minimizes correlations between differ- 
ent Wilson chains and in which, it turns out, the scat- 
tering matrix becomes diagonal. Moreover, the DMRG 
sweeping procedure allows the dimensions D of the MPS 
matrices in the system to be adjusted very flexibly. In- 
deed, in our case it was possible to reduce the matrix 
dimensions along the chain away from the dot consider- 
ably. The combination of all these resource-saving fea- 
tures makes it feasible to calculate the ground-state prop- 
erties of generic complex quantum impurity models using 
only relatively moderate numerical resources. 

The calculation of dynamical quantities like local spec- 
tral functions is, in principle, also possible for the star ge- 
ometry, for example, by suitably modifying the approach 
of Ref. to the present geometry. However, we expect 
that the increased computational costs of DMRG rela- 
tive to NRG for calculating dynamical quantities would 
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in this case likely offset the advantages of the star geom- 
etry. 

In closing, we would like to make the following com- 
ment: while we expect that a rotation to an optimal basis 
as described above should be applicable to a large class 
of impurity models, there may be cases where it does not 
work. In particular, we suspect that this might be the 
case for some models showing non-Fermi-liquid behav- 
ior, such as the two-channel spin-i Kondo model, where 
overscreening of the impurity is likely to lead to strong 
mutual correlations between all Wilson chains. A quan- 
titative analysis of this problem using the present star 
geometry approach is beyond the scope of the present in- 
vestigation but would be an interesting subject for future 



study. 
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